set.seed(42)
# dependencies
library(tidyverse)
library(metafor)
library(knitr)
library(kableExtra)
library(psych)
library(ggExtra)
library(scales)
# function to round all numeric vars in a data frame
round_df <- function(df, n_digits = 3) {
df %>% mutate_if(is.numeric, round, digits = n_digits)
}
probability_of_superiority_function <- function(x, y) {
nx <- length(x)
ny <- length(y)
rx <- sum(rank(c(x, y))[1:nx])
A = (rx / nx - (nx + 1) / 2) / ny
return(A)
}
spearman_brown_correction <- function(r_pearson) {
inversion <- ifelse(r_pearson < 0, -1, 1)
val <- round(2*abs(r_pearson)/(1 + abs(r_pearson)), 2)
return(val*inversion)
}
apa_p_value <- function(p){
p_formatted <- ifelse(p >= 0.001, paste("=", round(p, 3)),
ifelse(p < 0.001, "< .001", NA))
p_formatted <- gsub(pattern = "0.", replacement = ".", x = p_formatted, fixed = TRUE)
p_formatted
}
add_heterogeneity_metrics_to_forest <- function(fit) {
## more detailed
# bquote(paste("RE Model (Q(df = ", .(formatC(fit$k - 1)), ") = ",
# .(formatC(round(fit$QE, 2))),
# ", p ", .(formatC(apa_p_value(fit$QEp))),
# ", ", tau^2, " = ", .(formatC(round(fit$tau2, 1))),
# ", ", I^2, " = ", .(formatC(round(fit$I2, 1))),
# "%, ", H^2," = ", .(formatC(round(fit$H2, 1))), ")"))
## less detailed
bquote(paste("RE Model (", tau^2, " = ", .(formatC(round(fit$tau2, 1))),
", ", I^2, " = ", .(formatC(round(fit$I2, 1))),
"%, ", H^2," = ", .(formatC(round(fit$H2, 1))), ")"))
}
# table format in output
options(knitr.table.format = "html")
# create directory needed to save output
dir.create("models")
dir.create("plots")data_demographics <-
read_csv("../data/effect sizes for meta-analyses/data_demographics.csv")
data_test_retest_icc <-
read_csv("../data/effect sizes for meta-analyses/data_test_retest_icc.csv")
data_internal_consistency_oddeventrials <-
read_csv("../data/effect sizes for meta-analyses/data_internal_consistency_oddeventrials.csv")
data_internal_consistency_firstsecondhalves <-
read_csv("../data/effect sizes for meta-analyses/data_internal_consistency_firstsecondhalves.csv")
data_internal_consistency_permuted_estimates <-
read_csv("../data/effect sizes for meta-analyses/data_internal_consistency_permuted_estimates.csv") data_demographics %>%
summarize(mean_age = mean(age, na.rm = TRUE),
sd_age = sd(age, na.rm = TRUE)) %>%
round_df(1) %>%
kable() %>%
kable_styling(bootstrap_options = c("striped", "hover", "condensed"), full_width = FALSE)| mean_age | sd_age |
|---|---|
| 21 | 5.8 |
data_demographics %>%
count(gender) %>%
kable() %>%
kable_styling(bootstrap_options = c("striped", "hover", "condensed"), full_width = FALSE)| gender | n |
|---|---|
| Female | 318 |
| Male | 186 |
| Other | 1 |
| NA | 1145 |
data_test_retest_icc %>%
summarize(total_trt_n = sum(n)) %>%
kable() %>%
kable_styling(bootstrap_options = c("striped", "hover", "condensed"), full_width = FALSE)| total_trt_n |
|---|
| 340 |
data_internal_consistency_permuted_estimates %>%
summarize(total_ic_n = sum(n)) %>%
kable() %>%
kable_styling(bootstrap_options = c("striped", "hover", "condensed"), full_width = FALSE)| total_ic_n |
|---|
| 1546 |
Three methods of calculating internal consistency for the Implicit Relational Assessment Procedure (IRAP): - Odd vs even trials (mostly commonly reported in literature) - First vs second half of task (more common for the most popular implicit measure, the IAT; would allow a like-for-like comparison) - Cronbach’s alpha via permutation (most robust)
The IRAP effect is most often quantified using the D1 scoring algorithm (Greenwald, Banaji & Nosek, 2003). i.e., trim all RTs > 10,000 ms, then D1 = (mean RT in the inconsistent blocks - mean RT in the consistent blocks)/SD of RTs in both block types. D1 scores are used throughout. Because of this, it wasn’t possible to use Sam Parson’s splithalf R package. A workflow for permutation-based alpha estimation using D scores is provided below.
In each case, correlations between D scores calculated from the split halves are calculated, and then the Spearman-Brown prophecy formula is applied, i.e., r_sb = 2*r_pearson / (1+r_pearson). Importantly, this correction is applied to the absolute value of the pearson’s correlation, and then its sign is corrected to be congruent with the pearson’s correlation. This ensures that r_sb has a possible range of -1 to +1, as correlations should (whereas correction of native negative correlations allow for a lower limit of -Inf). I haven’t seen this treatment of SB corrections applied to negative correlations well explicated elsewhere, so it is important to note here.
For meta analysis, spearman brown correlations are converted to Fischers r-to-z for meta analysis. Estimates of the meta effect are then converted back for reporting. Heterogeneity metrics represent heterogeneity in Fischer’s r-to-z estimates. See here.
Bonett transformations
# fit random Effects model
fit_oddeven <- data_internal_consistency_oddeventrials %>%
rma(yi = yi,
vi = vi,
data = .,
slab = domain)
# make predictions
predictions_oddeven <-
predict(fit_oddeven, digits = 5) %>%
as.data.frame() %>%
round_df(2)
# plot
metafor::forest(fit_oddeven,
xlab = "Spearman-Brown correlation",
addcred = TRUE,
refline = FALSE,
transf = transf.iabt,
xlim = c(-1.4, 1.6),
at = c(0, 0.25, 0.5, 0.75, 1),
mlab = add_heterogeneity_metrics_to_forest(fit_oddeven))
text(-1.4, 39, "Study", pos = 4)
text(1.6, 39, "Spearman-Brown correlation [95% CI]", pos = 2)# summarize results
meta_effect <-
paste0("Meta analysis: k = ", fit_oddeven$k,
", r_sb = ", round(transf.iabt(predictions_oddeven$pred), 2),
", 95% CI [", round(transf.iabt(predictions_oddeven$ci.lb), 2), ", ",
round(transf.iabt(predictions_oddeven$ci.ub), 2), "]",
", 95% CR [", round(transf.iabt(predictions_oddeven$cr.lb), 2), ", ",
round(transf.iabt(predictions_oddeven$cr.ub), 2), "]")
meta_heterogeneity <-
paste0("Heterogeneity tests: Q(df = ", fit_oddeven$k - 1, ") = ", round(fit_oddeven$QE, 2),
", p ", ifelse(fit_oddeven$QEp < 0.0001, "< .0001", paste0("= ", as.character(round(fit_oddeven$QEp, 4)))),
", tau^2 = ", round(fit_oddeven$tau2, 2),
", I^2 = ", round(fit_oddeven$I2, 2),
"%, H^2 = ", round(fit_oddeven$H2, 2))Meta effect: Meta analysis: k = 37, r_sb = 0.54, 95% CI [0.47, 0.59], 95% CR [0.3, 0.7].
Heterogeneity: Heterogeneity tests: Q(df = 36) = 56.34, p = 0.0166, tau^2 = 0.04, I^2 = 29.04%, H^2 = 1.41.
Bonett transformations
# fit random Effects model
fit_firstsecond <- data_internal_consistency_firstsecondhalves %>%
rma(yi = yi,
vi = vi,
data = .,
slab = domain)
# make predictions
predictions_firstsecond <-
predict(fit_firstsecond, digits = 5) %>%
as.data.frame() %>%
round_df(2)
# plot
metafor::forest(fit_firstsecond,
xlab = "Spearman-Brown correlation",
addcred = TRUE,
refline = FALSE,
transf = transf.iabt,
xlim = c(-1.4, 1.6),
at = c(0, 0.25, 0.5, 0.75, 1),
mlab = add_heterogeneity_metrics_to_forest(fit_firstsecond))
text(-1.4, 39, "Study", pos = 4)
text(1.6, 39, "Spearman-Brown correlation [95% CI]", pos = 2)# summarize results
meta_effect <-
paste0("Meta analysis: k = ", fit_firstsecond$k,
", r_sb = ", round(transf.iabt(predictions_firstsecond$pred), 2),
", 95% CI [", round(transf.iabt(predictions_firstsecond$ci.lb), 2), ", ",
round(transf.iabt(predictions_firstsecond$ci.ub), 2), "]",
", 95% CR [", round(transf.iabt(predictions_firstsecond$cr.lb), 2), ", ",
round(transf.iabt(predictions_firstsecond$cr.ub), 2), "]")
meta_heterogeneity <-
paste0("Heterogeneity tests: Q(df = ", fit_firstsecond$k - 1, ") = ", round(fit_firstsecond$QE, 2),
", p ", ifelse(fit_firstsecond$QEp < 0.0001, "< .0001", paste0("= ", as.character(round(fit_firstsecond$QEp, 4)))),
", tau^2 = ", round(fit_firstsecond$tau2, 2),
", I^2 = ", round(fit_firstsecond$I2, 2),
"%, H^2 = ", round(fit_firstsecond$H2, 2))Meta effect: Meta analysis: k = 37, r_sb = 0.52, 95% CI [0.47, 0.57], 95% CR [0.47, 0.57].
Heterogeneity: Heterogeneity tests: Q(df = 36) = 40.09, p = 0.2938, tau^2 = 0, I^2 = 0%, H^2 = 1.
ie an estimate of alpha through large number number of random half splits (see Parsons, Kruijt, & Fox. 2019). Estimate of alpha obtained via mean of the resampled Spearman-Brown corrected split half correlations, when calculating D1 scores for each half and sampling 5000 permutations.
Bonett transformations
# fit random Effects model
fit_internal_consistency_permuted_estimates <- data_internal_consistency_permuted_estimates %>%
rma(yi = yi,
vi = vi,
data = .,
slab = domain)
# make predictions
predictions_permutations <-
predict(fit_internal_consistency_permuted_estimates, digits = 5) %>%
as.data.frame() %>%
round_df(2)
# plot
metafor::forest(fit_internal_consistency_permuted_estimates,
xlab = bquote(paste("Cronbach's ", alpha)),
addcred = TRUE,
refline = FALSE,
transf = transf.iabt,
xlim = c(-1.4, 1.6),
at = c(0, 0.25, 0.5, 0.75, 1),
mlab = add_heterogeneity_metrics_to_forest(fit_internal_consistency_permuted_estimates))
text(-1.4, 39, "Study", pos = 4)
text(1.6, 39, bquote(paste("Cronbach's ", alpha, " [95% CI]")), pos = 2)# summarize results
meta_effect <-
paste0("Meta analysis: k = ", fit_internal_consistency_permuted_estimates$k,
", alpha = ", round(transf.iabt(predictions_permutations$pred), 2),
", 95% CI [", round(transf.iabt(predictions_permutations$ci.lb), 2), ", ",
round(transf.iabt(predictions_permutations$ci.ub), 2), "]",
", 95% CR [", round(transf.iabt(predictions_permutations$cr.lb), 2), ", ",
round(transf.iabt(predictions_permutations$cr.ub), 2), "]")
meta_heterogeneity <-
paste0("Heterogeneity tests: Q(df = ", fit_internal_consistency_permuted_estimates$k - 1, ") = ", round(fit_internal_consistency_permuted_estimates$QE, 2),
", p ", ifelse(fit_internal_consistency_permuted_estimates$QEp < 0.0001, "< .0001", paste0("= ", as.character(round(fit_internal_consistency_permuted_estimates$QEp, 4)))),
", tau^2 = ", round(fit_internal_consistency_permuted_estimates$tau2, 2),
", I^2 = ", round(fit_internal_consistency_permuted_estimates$I2, 2),
"%, H^2 = ", round(fit_internal_consistency_permuted_estimates$H2, 2))
# save to disk for making pdf plot
write_rds(fit_internal_consistency_permuted_estimates, "models/fit_internal_consistency_permuted_estimates.rds")Meta effect: Meta analysis: k = 37, alpha = 0.52, 95% CI [0.47, 0.57], 95% CR [0.47, 0.57].
Heterogeneity: Heterogeneity tests: Q(df = 36) = 39.27, p = 0.3257, tau^2 = 0, I^2 = 0%, H^2 = 1.
Sexuality IRAP seems to be an outlier
if(!file.exists("models/fit_gosh_ic.rds")){
fit_gosh_ic <- gosh(fit_internal_consistency_permuted_estimates,
subsets = 5000,
parallel = "multicore", # change to "snow" on windows
ncpus = 4)
write_rds(fit_gosh_ic, "models/fit_gosh_ic.rds")
} else {
fit_gosh_ic <- read_rds("models/fit_gosh_ic.rds")
}
data_for_plotting_ic <-
bind_cols(as.tibble(fit_gosh_ic$incl) %>%
# the study to be separated is here listed by order (counting down from top of forest plot)
select(include = "18"),
as.tibble(fit_gosh_ic$res)) %>%
mutate(IRAP = ifelse(include, "Sexuality", "All other domains"),
alpha = transf.iabt(estimate)) %>%
# randomize rows to improve plotting
sample_frac(size = 1)
# ggplot(data_for_plotting_ic, aes(alpha, I2)) +
# geom_point(alpha = 0.4) +
# scale_color_viridis_d(begin = 0.3, end = 0.7, direction = -1) +
# theme_classic()
p1 <-
ggplot(data_for_plotting_ic, aes(alpha, I2, color = IRAP)) +
geom_point(alpha = 0.5) +
scale_color_viridis_d(begin = 0.35, end = 0.65, direction = -1) +
theme_classic() +
xlab(bquote("Cronbach's " ~alpha)) +
ylab(expression(italic("I")^"2")) +
theme(legend.position = c(0.2, 0.85),
legend.background = element_rect(fill = scales::alpha("white", 0)))
plot_gosh_ic <-
ggMarginal(p1,
type = "densigram",
size = 4,
margins = "both",
groupFill = TRUE)
plot_gosh_icwrite_rds(plot_gosh_ic, "models/plot_gosh_ic.rds")
data_for_plotting_ic %>%
group_by(IRAP) %>%
summarize(alpha_median = median(alpha, na.rm = TRUE),
alpha_ci_lower = quantile(alpha, 0.025, na.rm = TRUE),
alpha_ci_upper = quantile(alpha, 0.975, na.rm = TRUE),
I2_median = median(I2, na.rm = TRUE),
I2_ci_lower = quantile(I2, 0.025, na.rm = TRUE),
I2_ci_upper = quantile(I2, 0.975, na.rm = TRUE),
H2_median = median(H2, na.rm = TRUE),
H2_ci_lower = quantile(H2, 0.025, na.rm = TRUE),
H2_ci_upper = quantile(H2, 0.975, na.rm = TRUE)) %>%
round_df(2) %>%
kable() %>%
kable_styling(bootstrap_options = c("striped", "hover", "condensed"), full_width = FALSE)| IRAP | alpha_median | alpha_ci_lower | alpha_ci_upper | I2_median | I2_ci_lower | I2_ci_upper | H2_median | H2_ci_lower | H2_ci_upper |
|---|---|---|---|---|---|---|---|---|---|
| All other domains | 0.53 | 0.45 | 0.60 | 3.07 | 0.00 | 29.77 | 1.03 | 1.00 | 1.42 |
| Sexuality | 0.58 | 0.52 | 0.64 | 50.95 | 28.09 | 69.73 | 2.04 | 1.39 | 3.30 |
excluding Sexuality IRAP as an outlier
# fit random Effects model
fit_internal_consistency_permuted_estimates_sensitivity <- data_internal_consistency_permuted_estimates %>%
filter(!str_detect(domain, "Sexuality")) %>%
rma(yi = yi,
vi = vi,
data = .,
slab = domain)
# make predictions
predictions_permutations_sensitivity <-
predict(fit_internal_consistency_permuted_estimates_sensitivity, digits = 5) %>%
as.data.frame() %>%
round_df(2)
# plot
metafor::forest(fit_internal_consistency_permuted_estimates_sensitivity,
xlab = bquote(paste("Cronbach's ", alpha)),
addcred = TRUE,
refline = FALSE,
transf = transf.iabt,
xlim = c(-1.4, 1.6),
at = c(0, 0.25, 0.5, 0.75, 1),
mlab = add_heterogeneity_metrics_to_forest(fit_internal_consistency_permuted_estimates_sensitivity))
text(-1.4, 38, "Study", pos = 4)
text(1.6, 38, bquote(paste("Cronbach's ", alpha, " [95% CI]")), pos = 2)# summarize results
meta_effect <-
paste0("Meta analysis: k = ", fit_internal_consistency_permuted_estimates_sensitivity$k,
", alpha = ", round(transf.iabt(predictions_permutations_sensitivity$pred), 2),
", 95% CI [", round(transf.iabt(predictions_permutations_sensitivity$ci.lb), 2), ", ",
round(transf.iabt(predictions_permutations_sensitivity$ci.ub), 2), "]",
", 95% CR [", round(transf.iabt(predictions_permutations_sensitivity$cr.lb), 2), ", ",
round(transf.iabt(predictions_permutations_sensitivity$cr.ub), 2), "]")
meta_heterogeneity <-
paste0("Heterogeneity tests: Q(df = ", fit_internal_consistency_permuted_estimates_sensitivity$k - 1, ") = ", round(fit_internal_consistency_permuted_estimates_sensitivity$QE, 2),
", p ", ifelse(fit_internal_consistency_permuted_estimates_sensitivity$QEp < 0.0001, "< .0001", paste0("= ", as.character(round(fit_internal_consistency_permuted_estimates_sensitivity$QEp, 4)))),
", tau^2 = ", round(fit_internal_consistency_permuted_estimates_sensitivity$tau2, 2),
", I^2 = ", round(fit_internal_consistency_permuted_estimates_sensitivity$I2, 2),
"%, H^2 = ", round(fit_internal_consistency_permuted_estimates_sensitivity$H2, 2))
# save to disk for making pdf plot
write_rds(fit_internal_consistency_permuted_estimates_sensitivity, "models/fit_internal_consistency_permuted_estimates_sensitivity.rds")Meta effect: Meta analysis: k = 36, alpha = 0.51, 95% CI [0.46, 0.56], 95% CR [0.46, 0.56].
Heterogeneity: Heterogeneity tests: Q(df = 35) = 23.07, p = 0.9392, tau^2 = 0, I^2 = 0%, H^2 = 1.
Sexuality IRAP seems to be an outlier
if(!file.exists("models/fit_gosh_ic_sensitivity.rds")){
fit_gosh_ic_sensitivity <- gosh(fit_internal_consistency_permuted_estimates_sensitivity,
subsets = 5000,
parallel = "multicore", # change to "snow" on windows
ncpus = 4)
write_rds(fit_gosh_ic_sensitivity, "models/fit_gosh_ic_sensitivity.rds")
} else {
fit_gosh_ic_sensitivity <- read_rds("models/fit_gosh_ic_sensitivity.rds")
}
data_for_plotting_ic_sensitivity <-
as.tibble(fit_gosh_ic_sensitivity$res) %>%
mutate(Type = "Other IRAPs",
alpha = transf.iabt(estimate)) %>%
# randomize rows to improve plotting
sample_frac(size = 1)
# ggplot(data_for_plotting_ic_sensitivity, aes(alpha, I2)) +
# geom_point(alpha = 0.4) +
# scale_color_viridis_d(begin = 0.3, end = 0.7, direction = -1) +
# theme_classic()
p3 <-
ggplot(data_for_plotting_ic_sensitivity, aes(alpha, I2, color = Type)) +
geom_point(alpha = 0.5) +
scale_color_viridis_d(begin = 0.35, end = 0.65, direction = -1) +
theme_classic() +
xlab(bquote("Cronbach's " ~alpha)) +
ylab(expression(italic("I")^"2")) +
theme(legend.position = c(0.2, 0.85)) +
theme(legend.position = "none")
plot_gosh_ic_sensitivity <-
ggMarginal(p3,
type = "densigram",
size = 4,
margins = "both",
groupFill = TRUE)
plot_gosh_ic_sensitivitywrite_rds(plot_gosh_ic_sensitivity, "models/plot_gosh_ic_sensitivity.rds")
data_for_plotting_ic_sensitivity %>%
summarize(alpha_median = median(alpha, na.rm = TRUE),
alpha_ci_lower = quantile(alpha, 0.025, na.rm = TRUE),
alpha_ci_upper = quantile(alpha, 0.975, na.rm = TRUE),
I2_median = median(I2, na.rm = TRUE),
I2_ci_lower = quantile(I2, 0.025, na.rm = TRUE),
I2_ci_upper = quantile(I2, 0.975, na.rm = TRUE),
H2_median = median(H2, na.rm = TRUE),
H2_ci_lower = quantile(H2, 0.025, na.rm = TRUE),
H2_ci_upper = quantile(H2, 0.975, na.rm = TRUE)) %>%
round_df(2) %>%
kable() %>%
kable_styling(bootstrap_options = c("striped", "hover", "condensed"), full_width = FALSE)| alpha_median | alpha_ci_lower | alpha_ci_upper | I2_median | I2_ci_lower | I2_ci_upper | H2_median | H2_ci_lower | H2_ci_upper |
|---|---|---|---|---|---|---|---|---|
| 0.53 | 0.45 | 0.6 | 3.32 | 0 | 29.81 | 1.03 | 1 | 1.42 |
I.e., Intraclass Correlation Coefficients.
“It has been argued that test-retest reliability should reflect absolute agreement, rather than consistency, between measurements (Koo & Li, 2016). For example, a perfect correlation between scores at two time points may occur also when there is a systematic difference between time points (i.e. a difference that is about equal for all participants).” (Parsons, Kruijt, & Fox, 2019). Absolute agreement therefore also takes within-participant changes into account in its denominator, where consistency does not.
# fit random Effects model
fit_test_retest_icc <- data_test_retest_icc %>%
rma(yi = ICC,
sei = se,
data = .,
slab = domain)
# make predictions
predictions_test_retest_icc <-
predict(fit_test_retest_icc, digits = 5) %>%
as.data.frame() %>%
round_df(2)
# plot
metafor::forest(fit_test_retest_icc,
xlab = "Absolute Agreement ICC",
addcred = TRUE,
refline = FALSE,
xlim = c(-1.5, 1.7),
at = c(-0.5, -0.25, 0, 0.25, 0.5, 0.75, 1),
mlab = add_heterogeneity_metrics_to_forest(fit_test_retest_icc))
text(-1.5, 10, "Study", pos = 4)
text(1.7, 10, "ICC [95% CI]", pos = 2)# summarize results
meta_effect <-
paste0("Meta analysis: k = ", fit_test_retest_icc$k,
", ICC2 = ", round(transf.ztor(predictions_test_retest_icc$pred), 2),
", 95% CI [", round(transf.ztor(predictions_test_retest_icc$ci.lb), 2), ", ",
round(transf.ztor(predictions_test_retest_icc$ci.ub), 2), "]",
", 95% CR [", round(transf.ztor(predictions_test_retest_icc$cr.lb), 2), ", ",
round(transf.ztor(predictions_test_retest_icc$cr.ub), 2), "]")
meta_heterogeneity <-
paste0("Heterogeneity tests: Q(df = ", fit_test_retest_icc$k - 1, ") = ", round(fit_test_retest_icc$QE, 2),
", p ", ifelse(fit_test_retest_icc$QEp < 0.0001, "< .0001", paste0("= ", as.character(round(fit_test_retest_icc$QEp, 4)))),
", tau^2 = ", round(fit_test_retest_icc$tau2, 2),
", I^2 = ", round(fit_test_retest_icc$I2, 2),
"%, H^2 = ", round(fit_test_retest_icc$H2, 2))
# write to disk for pdf plots
write_rds(fit_test_retest_icc, "models/fit_test_retest_icc.rds")Meta effect: Meta analysis: k = 8, ICC2 = 0.18, 95% CI [0.05, 0.3], 95% CR [-0.11, 0.44].
Heterogeneity: Heterogeneity tests: Q(df = 7) = 15.89, p = 0.0261, tau^2 = 0.02, I^2 = 54.92%, H^2 = 2.22.
if(!file.exists("models/fit_gosh_trt.rds")){
fit_gosh_trt <- gosh(fit_test_retest_icc,
parallel = "multicore", # change to "snow" on windows
ncpus = 4)
write_rds(fit_gosh_trt, "models/fit_gosh_trt.rds")
} else {
fit_gosh_trt <- read_rds("models/fit_gosh_trt.rds")
}
data_for_plotting_trt <-
as.tibble(fit_gosh_trt$res) %>%
mutate(Type = "Other IRAPs",
ICC = transf.ztor(estimate)) %>%
# randomize rows to improve plotting
sample_frac(size = 1)
# ggplot(data_for_plotting_ic, aes(alpha, I2)) +
# geom_point(alpha = 0.4) +
# scale_color_viridis_d(begin = 0.3, end = 0.7, direction = -1) +
# theme_classic()
p2 <-
ggplot(data_for_plotting_trt, aes(ICC, I2, color = Type)) +
geom_point(alpha = 0.8) +
scale_color_viridis_d(begin = 0.35, end = 0.65, direction = -1) +
theme_classic() +
xlab("Absolute Agreement ICC") +
ylab(expression(italic("I")^"2")) +
theme(legend.position = "none")
plot_gosh_trt <-
ggMarginal(p2,
type = "densigram",
size = 4,
margins = "both",
groupFill = TRUE)
plot_gosh_trtwrite_rds(plot_gosh_trt, "models/plot_gosh_trt.rds")
data_for_plotting_trt %>%
summarize(ICC_median = median(ICC, na.rm = TRUE),
ICC_ci_lower = quantile(ICC, 0.025, na.rm = TRUE),
ICC_ci_upper = quantile(ICC, 0.975, na.rm = TRUE),
I2_median = median(I2, na.rm = TRUE),
I2_ci_lower = quantile(I2, 0.025, na.rm = TRUE),
I2_ci_upper = quantile(I2, 0.975, na.rm = TRUE),
H2_median = median(H2, na.rm = TRUE),
H2_ci_lower = quantile(H2, 0.025, na.rm = TRUE),
H2_ci_upper = quantile(H2, 0.975, na.rm = TRUE)) %>%
round_df(2) %>%
kable() %>%
kable_styling(bootstrap_options = c("striped", "hover", "condensed"), full_width = FALSE)| ICC_median | ICC_ci_lower | ICC_ci_upper | I2_median | I2_ci_lower | I2_ci_upper | H2_median | H2_ci_lower | H2_ci_upper |
|---|---|---|---|---|---|---|---|---|
| 0.18 | 0 | 0.32 | 53.34 | 0 | 75.74 | 2.14 | 1 | 4.12 |
Maximum observed correlation of two measures is a function of their true correlation and the reliability (i.e., self-correlation) of each measure.
\(r_{xy}=\frac{\rho_{xy}}{\sqrt{\rho_{xx}\rho_{yy}}}\)
observed_r <- function(true_r, reliability_a, reliability_b) {
round(true_r * sqrt(reliability_a*reliability_b), 2)
}
max_cor <- observed_r(true_r = 1.0, reliability_a = transf.iabt(predictions_permutations_sensitivity$pred), reliability_b = 1.0)
medium_cor <- observed_r(true_r = 0.5, reliability_a = transf.iabt(predictions_permutations_sensitivity$pred), reliability_b = 1.0)Given this meta estimate of reliability, assuming that the internal consistency of an external variable is perfect (1.0), the maximum correlation between the IRAP and that external variable (i.e., when true correlation is 1.0) is 0.71. For a true correlation of medium size (r = .5), the observed correlation would be 0.36.
max_cor <- observed_r(true_r = 1.0, reliability_a = transf.ztor(predictions_test_retest_icc$pred), reliability_b = 1.0)
medium_cor <- observed_r(true_r = 0.5, reliability_a = transf.iabt(predictions_test_retest_icc$pred), reliability_b = 1.0)Given this meta estimate of reliability, assuming that the internal consistency of an external variable is perfect (1.0), the maximum correlation between the IRAP and that external variable (i.e., when true correlation is 1.0) is 0.42. For a true correlation of medium size (r = .5), the observed correlation would be 0.2.
The Spearman-Brown prediction formula can be rearranged to make prediction about how the length of the test would need to change to meet a goal reliability:
\(\rho_{xx'}^*=\frac {n\rho_{xx'}}{1+(n-1)\rho_{xx'}}\)
where \(\rho_{xx'}^*\) refers to the goal reliability, \(\rho_{xx'}\) refers to the current reliability, and \(n\) refers to the test length multiplier.
current_ic <- transf.iabt(predictions_permutations_sensitivity$pred)
goal_ic <- 0.70
ic_length_increase_.70 <- round((goal_ic*(1 - current_ic)) / (current_ic*(1 - goal_ic)), 1)
goal_ic <- 0.80
ic_length_increase_.80 <- round((goal_ic*(1 - current_ic)) / (current_ic*(1 - goal_ic)), 1)
goal_ic <- 0.90
ic_length_increase_.90 <- round((goal_ic*(1 - current_ic)) / (current_ic*(1 - goal_ic)), 1)Using \(a\) = 0.51, the IRAP would have to be 2.3 times longer for it to have an internal consistency of >=.70, 3.9 times longer for it to have an internal consistency of >=.80, or 8.7 times longer for it to have an internal consistency of >=.90.
The Spearman-Brown prediction formula can be rearranged to make prediction about how the length of the test would need to change to meet a goal reliability:
\({\rho }_{xx'}^{*}={\frac {n{\rho }_{xx'}}{1+(n-1){\rho }_{xx'}}}\)
where \({\rho }_{xx'}^{*}\) refers to the goal reliability, \({\rho }_{xx'}\) refers to the current reliability, and \({n}\) refers to the test length multiplier.
current_test_retest <- transf.ztor(predictions_test_retest_icc$pred)
goal_test_retest <- 0.90
trt_length_increase_.90 <- round((goal_test_retest*(1 - current_test_retest)) / (current_test_retest*(1 - goal_test_retest)), 1)
goal_test_retest <- 0.80
trt_length_increase_.80 <- round((goal_test_retest*(1 - current_test_retest)) / (current_test_retest*(1 - goal_test_retest)), 1)
goal_test_retest <- 0.70
trt_length_increase_.70 <- round((goal_test_retest*(1 - current_test_retest)) / (current_test_retest*(1 - goal_test_retest)), 1)Using the absolute agreement estimate of 0.18, the IRAP would have to be 10.8 times longer for it to have a test-retest of >=.70, 18.5 times longer for it to have a test-retest of >=.80, or 41.5 times longer for it to have a test-retest of >=.90.